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CN Abstract 

I— I Plotting solution sets for particular equations may be complicated by the exis- 

"^] tence of turning points. Here we describe an algorithm which not only overcomes 

^^ such problematic points, but does so in the most general of settings. Applica- 

tions of the algorithm are highlighted through two examples: the first provides 
verification, while the second demonstrates a non-trivial application. The latter 
2 is followed by a thorough run-time analysis. While both examples deal with 

^ bivariate equations, it is discussed how the algorithm may be generalized for 

'""' space curves in R^. 
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C^ 1. Introduction 

q 

I- . In this paper we consider curves determined by equations of the form 



.fi^,y)^o (1) 



L| where /:/^M, /cMxMisa product of open intervals. Equations such as (1) 

• ^H are often referred to as implicit equations since in general there does not exist 

^^ an explicit, unique function g such that y = g{x). A canonical example is that 

'^ of the unit circle, with /(.t, y) = x^ + y^ — 1 = 0, which cannot be rearranged to 

isolate y as a function of x. The following discussion also applies to equations 
of the form /(x; a) = where a is a bifurcation parameter. 

The purpose of this paper is to address the problem of plotting solution 
curves to (1) with turning points (a precise definition will be presented shortly). 
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While a method for deahng with this problem has already been developed by 
Keller [7, 8] (via a pseudo-arc-length parametrization), it rehes on the turning 
points to be anything but cusps, and thus cannot be used on curves without 
some understanding of their profile. Conversely, the algorithm that we will 
present not only requires no prior knowledge of the solution curve, but also 
approaches the problem in what we consider to be a more intuitive manner. 

It should be said, however, that although we acknowledge that there are 
other methods to deal with this problem, we will make no attempt to compare 
them. Essentially, this paper is to be self-contained, and therefore our only 
concern is the explanation of our proposed algorithm, and its own benefits, not 
its relative benefits. 



2. Background 

Before a discussion of the algorithm may begin, we must first define the term 
turning point. 

Definition 1. Let / : M^ ^ M, and 5 = {{x,y) G M^ . f(x,y) = O}. A point 
{x*,y*) G S* is a turning point of f{x,y) = if {x*,y*) G S, and there exists 
S > such that at least one of the following statements is true: 

0<x-x* <S, 0<\y-y*\<S}nS = 0. 
0< |a;-a;*| <5, < y - y* < 6} n S ^ 0. 
0<x* ~x<d, 0<\y-y*\<S}r\S = 0. 
0<\x~x*\<5, 0<y* -y <S}r\S = 0. 
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Figure 1: Illustration of the different types of turning points 



Definition I may be envisioned in the following way: if one were to restrict 
themselves to a small enough neighborhood about {x*,y*), one of the open 
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Figure 2: Demonstration of algorithm on a Type 1 turning point. 

"half-balls" centered at {x*,y*) contained in this neighborhood would have no 
intersection with S (see Figure 1). 

The remainder of this paper will deal with the discrete set: 



'fe = {{(^..y.)},1i C K2 : Mx„y,) = 0, 1 < j < 7v} 



Sk = 



(2) 



where Sk is the approximation of S, and /j. is the discrete approximation /. For 
a detailed discussion of approximating turning points in discrete spaces see the 
review by Cliffe, Spence, and Tavener[9]. 



3. Algorithm 

Supposing that {xj,yj) G Sk has been identified as a turning point of any 
type, let R be the closed disk of radius r centered at {xj,yj). The algorithm is 
then performed in the following steps: 

Step 1. Uniformly scan the boundary of R, dR, for solutions to f{x,y) ~ 0. 

Denote the set of solutions as TZ. 
Step 2. (a) liTZ = 0, we may stop here: no other points in M^ satisfy f{x, y) = 
0. 
(b) liTZ^ 0, select {xi,yi) £ Sk, i < j, to be a "reference point" and 
select the minimal point (s', t') e TZ such that 0(s', t') < ip{s, t) for 
all (s,i) S TZ, where 



0(s,t) = 



1 



Step 3. (a) Create the vector v = (s' — Xj,t' — yj), 



(3) 



(b) Determine the largest axial component oiv, d — max | |s' — Xj|, |t' — j/j||. 

(c) Use d to determine the new direction of iteration. 

To clear up any possible ambiguity, we will elaborate on how one might 
implement the above strategy on a Type 1 turning point. 

Step 1. Since {xj,yj) is a Type 1 turning point, we know that there exists an 
open ball of radius 5 centered at {xj,yj) whose eastern hemisphere is disjoint 
from Sk- Therefore, if any solution to the equation f(x,y) — exists, it must 
exist in the western hemisphere. Hence, we describe 9i?+, our reduced search 
path, by 

dR+ = {{x,y) eM."^ -.{x- x^f + {y - yjf == (5^, a; < x,) 

which is easily seen to be parametrized by 

g{e) = {5 cos{e) +x J, 5 s\TL{e)+yj), e e [7r/2,37r/2] 

Thus, if n points are used to uniformly sample dR^, the mesh would look like 

g^ = (6cos( — {n + 2i)j + Xj,6sm( — {n + 2i)j + yA 

for i ~ 0,1, . . . ,n — 1. 

The points gi are then checked to see if f{gi) = 0. If gi is solution to 
f{x,y) = 0, it should be recorded in the set TZ. It should be noted that the 
method employed to check whether gi is a solution to f{x,y) = is irrelevant. 

Step 2. We have found that it is best to pick an {xi, yi) that is contained in R. 
In this way, when we minimize the cost function (j), see (3), we are searching for 
the the point in TZ which is furthest from our the reference point. If we assume 
{xj,yj) is not the final point of Sk, then the furthest point from the reference 
point is not expected to be any of the previous points in Sk ■ Picking a reference 
point that is "too far" from {xj,yj) can potentially lead to choosing {s',t') € TZ 
such that {s',t') is a previous solution. When the algorithm is implemented 
inside an iterated plotting program (one would do this for autonomous plotting), 
this has the consequence of leading the plotting program to retrace its old steps. 
From here it is only a matter of evaluating (j>(s,t) for all {s,t) £ TZ, and 
finding the smallest value. The coordinate which minimizes (j) shall be called 
{s',t'). While there is no garentee that there exists a single coordinate, such 
that (j) is minimized uniquely; in practice, we seldom arrive at two, or more 
points, that equally minimize (j). 

Step 3. Make the direction vector, v, from {xj,yj) to (s',t'). Then identify 
which component of v is larger in magnitude. If the y-componcnt is larger, 
begin iterating along the y-axis (the direction is determined by the sign of the 
y-component). 

Iteration should thence begin at {s',t'). This is to compensate for the fact 
that it may happen that the turning point was approached via iteration along 




Figure 3: Astroid plot generated using the turning point algorithm. Inset: parametrized curve 
(solid), approximated curve (circles) 



the positive x-direction, yet the algorithm has concluded that the new direction 
is along the negative x-direction. If we do not begin at (s', t'), the program will 
just back-track along already found points of Sk- 
it should be noted that while the algorithm has been presented in the two- 
dimensional case, it is easily generalized to apply to equations on M". For 
instance, if we were to consider the case /(x, y, z) = 0, then a simple modification 
of Definition 1 including half-spheres, instead of half-balls, would allow for this 
algorithm to work with space curves. 

4. Verification 

To demonstrate the effectiveness of the turning point algorithm, we chose to 
implement it on the astroid, a well known implicit equation given by f{x,y) = 



„2/3 



y 



^' "^ — 1 = 0. The resulting curve can be found in Figure 3. 



The accuracy of the curve generated with our algorithm was measured in 
percent error (4), and depended on several different parameters: r, the radius of 
dR+, k how many "steps" back the reference point was from the turning point, 
and n the grid size of 9i?+. Each parameter was varied until notable values 
were found. 



PercentError(xj , yj) 



y] 



-^ 



^/x'^ + : 



y 



X 100% 



(4) 



It was found that a radius r € [0.0001(5, 1005] (where 6 is the magnitude of 
the step-size A^ or Aj,), 1 < A: < 10 and 4 < n < 10 yielded a maximum percent 
error ^ 0.1%. However, for n < 4 the algorithm failed to correctly navigate the 
turning point, and the plot ceased to continue. 

5. Example: Lubrication Model 

Consider the problem of determining the thickness of a thin film of lubri- 
cant inside a rotating cylinder. Using the lubrication approximation to the 
Navier-Stokes equations developed in [10, 11], one may describe the steady- 
state thickness, h(9), by 

|(V + 0-ic„.,<,).f-l (5) 

where Q is the non-dimensional fiux and e is dependent on the surface tension 
and rate of rotation. Coupling (5) with 

M= / h{e)de (6) 

it is not unreasonable to ask how the curve of points satisfying (5) and (6) 
behaves for a particular e. Because this curve may be characterized by a solution 
set to some equation f^iQ, M) = 0, we will let f^ denote the curve. 

Since our intention is to demonstrate the robustness of the algorithm, and 
not the challenge behind solving fe{Q, M) = when either Q or M is fixed, we 
refer the curious reader to the works of Ashmore et at [14] and Benilov et al. 
[12] for detailed exposition of that task. 

Briefly, we discretized ft, on a uniform grid, and use Fourier spectral dif- 
ferentiation matrices which account for the 27r-periodicity of h. A MATLAB 
program was written to accept very small values of e and perform an iterated 
Newton-Raphson procedure to generate fe given a particular parameter and 
direction. Once a turning point was identified, the algorithm was performed. 

It was found that for e < 10^'^ complicated loop behavior similar to that 
exhibited in Figure 4 appeared. With the algorithm deployed to automate plot- 
ting, each successive turning point was over come with no difficulty. Interested 
readers may also refer to [12, 13, 14] for recent literature on this topic, as well 
as [15, 16, 17]. 

6. Run-Time Analysis 

Due to the various settings in which the algorithm may be applied, a general 
run-time analysis is impossible. Instead, we will analyze the case when / is not 
explicitly given, but inferred from differential or integral equations. In particu- 
lar, we will perform a run-time analysis on the Lubrication Model example. As 
demonstrated in Section 5, it is clear that these settings are very useful in the 
study of certain ordinary and partial differential equations. 
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Figure 4: A portion of the (Q, M)-bifurcation diagram generated for e ■ 
[18]. 
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In this particular instance, we had two particular equations (5) and (6) which 
allowed us to solve for Q and M given a function h : [0, 27r] — > M+. By approxi- 
mating h with hm (an to x 1 column vector), and taking the appropriate Taylor 
expansion of elements of (5), we were able to perform parameter continuation 
along the M-axis with a to x to. matrix approximation of the operator /. 

With similar linearization techniques we were able to perform parameter 
continuation along the Q-axis with a (to- + 1) x (m + 1) matrix approximation 
of /. (For more details about these techniques see Bcnilov [12].) 

In this setting, performing the Newton-Raphson method prescribed by Ash- 
more [14] and/or Benilov [12], and thus completing Step 1 in the algorithm, is 
bottle-necked by the speed of matrix multiplication. 

Now, an algorithm for matrix multiplication of ?Ti x to matrices has been 
found with run-time ©(tti^"^^^) (i.e. the Coppersmith- Winograd algorithm[19]). 
However, this algorithm is particularly hard to implement, and is usually used 
for theoretic bounds above all else. Thus, in the spirit of creating a practical 
run-time argument, we assume that Strassen's algorithm is used, and therefore 
the calculations in Step 1 can be carried out in ©(to^'*^^) [20]. Also, because 
Newton's method is usually performed in a for-loop for a fixed number of itera- 
tions, £ say. Step 1 may be preformed in 0{nllm'^-^^'^) ~ 0{nm?-^'^'^)^ where n 
is the number of sample points on dR+ . 

Now, because calculating (p is reasonably assumed to take constant time, it 
is clear that the most cumbersome part of Step 2 is determining the point which 
minimizes (j) (which is a sorting problem). Thus, if we assume TZ has n points 
in it (the maximal amount), then Step 2 may be performed in Q{n\ogn). This 
follows from the fact that [20] shows that that any comparison sort algorithm 



requires f2(nlogn) in the worst case. Hence, taking any version of heap or 
merge-sort to perform Step 2 will require Q{n\ogn). 

While there are a myriad of ways one may perform Step 3, any competent 
implementation should run in constant time since this step has only manipulates 
two points in K^. 

Hence, the algorithm can be found to run in C'(n?7i^®''^ + nlogn). However, 
as we saw in the astroid example, n need not exceed 10, or so, points. Thus, 
if we approximate h with any reasonable amount of resolution (i.e. ?n = 2^ or 
m = 2^), then n -^ m, and we may consider it something of a constant. Hence, 
the algorithm can be expected to run in ©(m^-^^^). 

7. Conclusion 

It is easy to see that turning points make automated plotting a challenge, 
if not impossible. To overcome this issue we have outlined an algorithm that 
determines the direction iteration should continue in after reaching a turning 
point. As demonstrated with the astroid and thin-film differential equation 
examples, the algorithm overcomes turning points accurately, and effectively, 
while asking little in terms of computational power. It should be said that 
implicit functions and dynamical systems far from exhaust the areas where 
this algorithm may be helpful. With the amount of research being done today 
throughout the various sciences (pure and applied alike), it is clear that this 
algorithm can be applied anywhere data visualization is possible. 
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